#!/usr/bin/python
#Shreads contigs file into fragments of desired sizes. 
#Example: use these reads as 'fake' reads in 454 assemblies.
#Usage: shredOverlappingReads.py seq.fasta seq.qual 400 100
#Needs sequence fasta to add quality values to it. 
#Number of fasta records should match number of qual records

import sys
import re
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

seqfile = sys.argv[1]
sf = SeqIO.parse(seqfile, "fasta")

qualfile = sys.argv[2]
qf = SeqIO.parse(qualfile, "qual")

size = int(sys.argv[3])
overlap = int(sys.argv[4])

start = 0
stop = size

for s, q in zip(sf, qf):
    start = 0
    stop = size
    while stop <= len(s.seq)+size:  
        sid = s.id + "_" + str(start) + "_" + str(start + len(s.seq[start:stop]))
        rec = q.letter_annotations["phred_quality"]
        record = SeqRecord(s.seq, id=sid, description="")
        record.letter_annotations["phred_quality"] = rec
        print record[start:stop].format("qual").strip("\n")
        if stop > len(s.seq):
            break #stop if last one is already printed
        start += size - overlap
        stop += size - overlap

